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Starting with the early alchemists, a holy grail of science has been to make desired materials by modifying the 
attributes of basic building blocks. Building blocks that show promise for assembling new complex materials 
can be synthesized at the nanoscale with attributes that would astonish the ancient alchemists in their versatility. 
However, this versatility means that making direct connection between building block attributes and bulk behavior 
is both necessary for rationally engineering materials, and difficult because building block attributes can be altered 
in many ways. Here we show how to exploit the malleability of the valence of colloidal nanoparticle “elements” 
to directly and quantitatively link building block attributes to bulk behavior through a statistical thermodynamic 
framework we term “digital alchemy”. We use this framework to optimize building blocks for a given target 
structure, and to determine which building block attributes are most important to control for self assembly, 
through a set of novel thermodynamic response functions, moduli and susceptibilities. We thereby establish direct 
links between the attributes of colloidal building blocks and the bulk structures they form. Moreover, our results 
give concrete solutions to the more general conceptual challenge of optimizing emergent behaviors in nature, and 
can be applied to other types of matter. As examples, we apply digital alchemy to systems of truncated tetrahedra, 
rhombic dodecahedra, and isotropically interacting spheres that self assemble diamond, FCC, and icosahedral 
quasicrystal structures, 


Mendeleev’s tabular organization of the elements (H 2 ] by 
atomic valence 01 has served for more than 140 years as 
a heuristic that relates properties of the atomic elements to 
how they arrange in bulk structures. However, attempts to 
understand how properties of bulk structures relate to atomic 
properties predate Mendeleev and, in fact, modern science (4), 
and are complicated by the fact that the chemical manipulation 
of atoms is prohibited by the quantization of both electrical 
charge and angular momentum. Fortunately for Mendeleev, 
this quantization constrains Nature to only about 80 stable 
elements, and limits elemental properties and bulk behaviors so 
that the elements can be tabulated by valence. In fact, starting 
with technetium 0 in the 1930s, new atomic elements have 
only been produced artificially (as suggested by the etymology 
of the name “technetium” 0) by a-particle bombardment, 
fusion, or other nuclear techniques that finally realized the 
ancient alchemists’ goal of transmuting the elements. 

In contrast, an inexhaustible array of new “elements” can be 
synthesized as patchy particles. mm However, the exploding 
diversity of patchy particles lISl-HOll or, more generally, col¬ 
loidal “elements” means that there are now so many types to 
synthesize and study that synthesizing them all and determin¬ 
ing their bulk behavior is no longer possible in practice. This 
fundamental impracticality means that, for colloid science to 
progress, scientists must first ask and answer the basic but 
daunting question “What elements should I make?” Materi¬ 
als science that starts with this question must be carried out 
in a fundamentally different way than traditional approaches, 
guided by the question: What is the optimal building block to 
make for a given structure, and why is it optimal? 

Constructing a periodic table of colloidal elements is easier 
said than done, however, because, unlike for atoms, colloid 
valence EHia is not discrete. Moreover, entropic colloid 
valence in na is a collective effect on that emerges only 


when colloids are crowded fl4l[T6l . A first step in constructing 
a periodic table of colloidal elements was taken by heuristically 
classifying building blocks according to their valence along 
anisotropy dimensions I® [14] that systematically and orthogo¬ 
nally vary colloidal element attributes. This sort of colloidal 
alchemy is now possible. 

Here we present a statistical thermodynamic framework that 
forms the basis for a new computational approach to building 
block design, which we term digital alchemy. Using this frame¬ 
work: (i) We show how to treat anisotropy dimensions HE) 
or other particle interaction parameters as thermodynamic vari¬ 
ables, and interpret their conjugate quantities. Treating particle 
interaction parameters thermodynamically means that the at¬ 
tributes of the colloidal “elements” we study can change so we 
refer to our methods as “alchemy” The term alchemy has been 
used previously in the modern era in the context of materials 
design, and these uses are either different in spirit from the 
present work ]H~8l , or are focused on computing global free 
energy differences in systems |[T9l 20]] in which intermediate 
state points are unphysical. A related investigation was also 
carried out in Ref. ED, which considered the effects of non- 
rigid colloid shape on crystallization mechanically, whereas 
here we study rigid colloids that fluctuate thermally. Though 
there are many systematic investigations of how particle shape 
or interactions affect structure lfT0l[T4ll22l[43lL we are aware 
of no work that attempts to directly probe the thermodynamic 
response of a system to a change in the attributes of its con¬ 
stituent building blocks, in analogy with pre-scientific attempts 
to modify chemical elements. fD (ii) We show how constitutive 
relations between anisotropy parameters and the thermody¬ 
namically conjugate variables we term “alchemical potentials” 
encode a broad class of detailed quantitative relations between 
building block attributes and bulk behavior. Further, we define 
new moduli and susceptibilites that describe stress-strain rela- 
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tionships between bulk structure and particle attributes, (iii) We 
show that these building block vs. bulk relationships persist 
in systems with entropy-driven, emergent collective behavior, 
(iv) We show how building block vs. bulk relationships can be 
used both to determine optimal particle shapes or interactions 
for given structures, and to compute the relative importance of 
different particle attributes for bulk behavior, (v) We report a 
detailed, general microscopic design rule for a macroscopic, 
entropy-driven, emergent behavior, (vi) We demonstrate this 
design rule in simulations that allow particle shape to fluctuate 
dynamically by showing that when particles are constrained to 
sit on a target lattice, they spontaneously adopt their preferred 
shape; that is, the shape that minimizes the free energy of the 
target structure at a given state point. 

Through all of these findings we demonstrate what the out¬ 
lines of a periodic table of colloidal elements might look like. 
In particular, because colloidal valence is not constrained, a 
colloidal periodic table cannot be as succinct as the atomic 
periodic table. However, because colloidal valence can be 
manipulated, a colloidal periodic table can encode detailed, 
quantitative relationships between building block attributes and 
bulk behavior, and can tell us what building blocks are optimal 
for a given structure, and why they are optimal. 

I. THEORETICAL RESULTS 

We consider a family of colloidal elements that can be 
described by a set of isotropic interaction potentials, or 
by anisotropy dimensions for enthalpic is or entropic m 
patches, with parameters {c^}. The particles are described by 
a classical Hamiltonian H that depends on the oti via a pair 
interaction between particles, and the rotational kinetic term in 
the Hamiltonian 


can be found with various methods. For brevity we start with 
the Shannon/Jaynes [46; 471 entropy 

s =-'Z[ 7Ta Mtiv) - /3(n a H - ( E ) 
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where we have set /^b = 1,7r a is the probability of finding the 
system in a state labelled a, (3 and /i* are Lagrange multipliers 
enforcing the thermal averages, N is the number of particles 
in the system (the factor of N is included here so that both 
and oti can be intensive quantities), and the summation should 
be interpreted schematically. Unless otherwise noted we will 
work in units where the particle volume = 1. To determine 
the partition function we maximize Eq. © with respect to 7r a . 
This gives, up to some normalization constant Z , 

^ (3) 

and fixing the normalization n a = 1 gives 

Z = ^2 ^ Nai ) . (4) 

<7 

We see that (3 = 1/T, the usual inverse temperature, and pi are 
generalized chemical potentials conjugate to the “charges” oti 
na that determine the building block attributes. To distinguish 
fii from the ordinary chemical potential, and since they act as 
sources for changing the “elemental” building blocks of the 
system, we refer to them as “alchemical” potentials. We define 
the thermodynamic potential for this ensemble as Z = e~ 
which gives 
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where p are momenta, L are angular momenta, / is the moment 
of inertia tensor, and U is the interaction potential that depends 
on particle positions q and orientations Q , and where we have 
suppressed particle indices. We consider systems in which the 
generalized particle coordinates and their conjugate momenta 
do not have explicit dependence on the oti. In this case, the 
oti have vanishing Poisson brackets with the Hamiltonian, and 
are invariants of the system: {c^, H} = 0. This is the case 
if, e.g ., a particle’s shape is independent of its generalized 
momentum and position. This would not hold, e.g., for systems 
with chemical gradients that cause a particle to swell in some 
locations more than others. Furthermore, we consider systems 
in which the ol{ themselves are mutually commuting, i.e., the 
order in which operations are applied to modify the building 
blocks is not important. We regard the a* as a set of mutually 
conserved charges, and it has been shown ll44|j45ll that there is a 
well-defined thermodynamic ensemble for any set of mutually 
commuting conserved charges. 

Formally, we consider a system where the fluctuate ther¬ 
mally about some averages (c^), and the energy fluctuates 
about an average (E). The partition function for this ensemble 


where ij is the packing fraction or density. This computes how 
the system responds to a change in alchemical potential, and 
in the thermodynamic limit (hereafter we will be concerned 
about the thermodynamic limit so we will drop the () notation) 
establishes a constitutive relation ai(rp T, {pj}). It is conve¬ 
nient to make a Legendre transformation F = (j> + JL piNai, 
and compute the constitutive relation gLi(iq, (3, {aj}) using the 
expression 


_ 1^ fdF\ 
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For notational simplicity, especially in cases where we con¬ 
sider a single c^, it will sometimes be convenient to drop the 
subscripts on a and fi. 

The constitutive relation fi{a) quantifies the thermodynamic 
response of a system to a change in the attributes of the con¬ 
stituent particles. See appendix for a discussion of higher 
order thermodynamic response functions. If the alchemical 
potential fi > 0 at some state point NVTa, then an infinites¬ 
imal increase in the alchemical parameter a would increase 
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the free energy of the system. Conversely if /x < 0 then an 
infinitesimal increase in a would decrease the free energy of 
the system. This has two important implications, (i) Locally 
optimal particle attributes a * are determined by the roots of 
the constitutive relation fi(a*) = 0 with positive slope. We 
show in appendix that the locations of these roots are invariant 
under reparametrizations of a. (ii) In hard particle systems, 
where the free energy simply measures the system entropy, 
/x directly measures how the number of states available to a 
system changes as a function of the particle shape, and so it 
can be used to systematically determine which particle features 
are most likely to come into contact, and provides explicit 
quantitative guidance on how to design shapes for structures. 
We demonstrate both of these implications below. 

In the next section, we explicitly compute /x in three ex¬ 
ample systems, and interpret the meaning and implications 
of each computation. We compute the constitutive relation 
/jLi(rj,T, {ay}) at 77 , T, {ay} numerically using Eq. ([ 6 ]) with 
the Bennett acceptance ratio method (48). Using this method, 
we compute /x at some {ay} by equilibrating several indepen¬ 
dent samples at nearby values ay + vhj, where v are constants 
chosen for an appropriate finite differencing scheme, and hj 
are finite differences. For a full description of the computation, 
see appendix. To determine valence for anisotropic particles, 
we use the potential of mean force and torque (PMFT), as 
described in Ref. m. 

In addition to constitutive relations between thermodynamic 
quantities (i.e. first order derivatives of the free energy), phys¬ 
ical systems are also frequently characterized by higher free 
energy derivatives: susceptibilities and moduli (see, e.g ., Refs. 
E3I2D). We define the alchemical modulus M a and suscepti¬ 
bility as 


M a 
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The extension to systems with several alchemical parameters 
is straightforward. We note that, like standard moduli (e.g. 
bulk, shear, Young’s), M a is a stress-strain relationship (50) . 
but the strain is in alchemical space rather than real space. 
Accordingly, alchemical modulus M a Eq. 0 at a* measures 
how sensitive the system is to deviations from the ideal particle 
properties. Similarly, like standard susceptibilities (e.g. com¬ 
pressibility) (50) . is a strain-stress relationship. Physically, 
e.g., by the fluctuation-dissipation theorem (see, e.g., Ref. (49)) 
Xa determines how quickly a system of, say, fluctuating shape 
relaxes when particles are perturbed from their equilibrium 
attributes. 


II. NUMERICAL RESULTS AND DISCUSSION 


We use our digital alchemy methodology to optimize build¬ 
ing blocks for self-assembly in three different case studies. The 
first two involve entropy-driven systems, which are among the 
most conceptually difficult in which to connect macroscopic 
and microscopic system properties because the macroscopic 
behaviors are intrinsically collective.! 14j f6i f5TH53) In the 


third study, we investigate an oscillating pair potential, which 
was recently shown [42| to self-assemble a one-component 
icosahedral quasicrystal, one of the most complex crystal struc¬ 
tures known. In each case, the details of the specific model 
are included in the discussion below. Details and extended 
discussion of the methods used in each case may be found in 
the SI. 


A. Truncated Tetrahedra 

We simulated a one-parameter family of truncated tetrahe¬ 
dra at moderate truncations known to self-assemble diamond 
lattices ED. We parametrized the truncation between a = 0 
(a tetrahedron maximally truncated so that it is an octahe¬ 
dron) and 1 (an untruncated, regular tetrahedron). With this 
parametrization (we discuss reparameterization invariance of 
our results in appendix) particles self-assembled diamond at 
a packing fraction of 77 = 0.6 between truncations of 0.25, 
and 0.475 (see Fig. |Tk). For reference, the Archimedean trun¬ 
cated tetrahedron (33) has a truncation of |. We performed 
standard Monte Carlo simulations (e.g. Ref. USD of systems 
of TV = 216 and 1000 particles at fixed volume. Polyhedra 
overlaps were checked using the GJK algorithm (55) . 

For the truncated tetrahedra, we computed the constitutive 
relation between vertex truncation a and its conjugate alchem¬ 
ical potential /x. We first computed fi in small systems of 
N = 216 particles, and found preliminary evidence for van¬ 
ishing alchemical potential (here, a free energy minimum) for 
0.35 < a* < 0.4, (Fig. [TJl, squares). To obtain higher preci¬ 
sion, and to test for finite size effects, we simulated systems 
of N = 1000 particles in the region surrounding the putative 
free energy minimum (Fig. [TJl, circles). From these alchemical 
potential computations we extracted the free energy of the sys¬ 
tem as a function of shape in the vicinity of the minimum (Fig. 
[ljl, inset), which we estimated by performing a weighted least 
squares fit to 


f3fi = (3M a (a — a*) , (8) 

from which we find the free energy minimum is at 

a* = 0.3736 ± 0.0001 , (9) 

and the alchemical modulus M a is 

/ 3M a (a = a*, 77 = 0.6) = 52.0 ± 0.3 . (10) 

We also constructed diamond densest packings (Fig. [TjO for 
truncated tetrahedra for all truncations (in increments of 0 . 001 ) 
at which self assembly into diamond lattices was reported in 
Ref. 133J, and find the curve has a maximum consistent with 
the Archimedean truncated tetrahedron at a a = 1/3. 

To directly examine the effects of shape modification on 
emergent valence USED, we computed the PMFT for sys¬ 
tems of N = 1000 truncated tetrahedra. For details of this 
computation, see Ref. 02). We computed the PMFT at a den¬ 
sity of 77 = 0.6 for a truncation of a = 0.25 (Fig. 0). and a 
truncation of a * (Fig. [ 2 J 5 ). The results for the first neighbor 
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Figure 1. Truncated tetrahedra at a range of truncations a (a) self-assemble a diamond lattice (b) <331 . A search for maximal diamond packing 
density, rid (f) would suggest optimal assembly at a a = 1/3, the Archimedean truncated tetrahedron. We compute the constitutive relation 
(d) ii(a) (Eq. ([ 6 ])) for hard truncated tetrahedra at 77 = 0.6. Squares (□) are results for systems with 216 particles, circles (O) for systems 
with 1000 particles; where visible, error bars are one standard deviation. The alchemical potential vanishes when the truncation is optimal for 
self-assembling diamond at this density when the truncation is approximately a* « 0.37. In the inset plot we reconstruct the free energy curve 
in the vicinity of the minimum. The increase in anisotropy a* above the geometric prediction a a arises because particles need to increase 
anisotropy to preserve tetrahedral valence at lower packing fractions, but if particles are too anisotropic, the simultaneous coordination of 
neighboring particles is sterically prohibited (c, see also Fig.|2|. We demonstrate this design rule (e; see also appendix Movie) by simulating 
tetrahedra with fluctuating shape at 11 — 0 in an Einstein crystal with spring constant k at packing fraction 77 = 0.6 and allowing the particles to 
find their optimal shape. The plot (e) shows that at low k the average truncation (a) is consistent with a* (N = 216 squares ■; N = 1000 
circles O). 


shell show particles have stronger tetrahedral valence at a* 
than at a = 0.25, which originates from the relatively larger 
hexagonal faces acting as stronger entropic patches fl4l . How¬ 
ever, we see that at a fluid density of 77 = 0.5, in the second 
neighbor shell when particles have the optimal truncation a * 
(cyan spots in Fig. [2]c), the next-to-nearest neighbors sit in an 
alternating arrangement, whereas the next-to-nearest neighbors 
for perfect tetrahedra (blue spots in Fig. [2]i ) are rotated by 
7 r/6. This indicates a non-alternating arrangement that coin¬ 
cides with polytetrahedral motifs not commensurate with the 


diamond lattice, which arises from steric constraints depicted 
in Fig. |T]c. To directly confirm this result, we performed sim¬ 
ulations in an NVTfi ensemble (i.e. both thermostated and 
“alchemostated”) to determine a(/i) at fi — 0 for N = 216 and 
1000 truncated tetrahedra with fluctuating shape in a diamond 
Einstein crystal. We initialized the system at low packing 
fraction 77 = 0.2 with fully truncated (i.e. octahedral, a = 0 ) 
particles, and slowly compressed the system to the target pack¬ 
ing fraction of 77 = 0 . 6 , after which we relaxed the spring 
constant. We observed that the process drove the particles to 
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Figure 2. Emergent valence encoded in the PMFT for truncated 
tetrahedra for a crystal at density 77 = 0.6 (a: a m 0.25, b: a = 
a * « 0.37), and a fluid at density 77 = 0.5 (c: a = a* « 0.37, 
d: a = 1.0). In the crystal we see that the particle at the optimal 
truncation a * (b) shows greater specificity of tetrahedral valence than 
at lower a (a), as expected. However, at fluid densities, we see that if 
the particle is too tetrahedral (d), the second neighbor shell is by 7r/6 
compared with lower truncations (c), and is incommensurate with the 
diamond lattice. 


spontaneously adopt a truncation consistent with our alchem¬ 
ical potential calculations at fixed shape. See Fig. [T£, and 




Figure 3. Rhombic dodecahedra have both four-fold and three-fold 
vertices (a). We determine the relative sensitivity to the truncation 
of each type of vertex by computing the constitutive relation 7 /( 77 ) 
(b) according to the (exaggerated) truncations shown in (a). We plot 
alchemical potentials for 4-fold truncations ( 7 / 4 , □ squares) and 3- 
fold truncations ( 7 / 3 , a triangles) at various densities for which the 
system self-assembles an fee lattice. We observed 7/4 < 773 < 0 at 
all densities, indicating that both vertex truncations improve assembly 
of the target crystal, but four-fold vertex truncations provide greater 
improvement. 


appendix movie. 

Our computation of the constitutive relation fi(a) for trun¬ 
cated tetrahedra that form a diamond lattice reveals several 
findings, (i) By determining that fi{a) has a root at a* ^0.37 
we have demonstrated that it is possible to find a thermodynam¬ 
ically optimal shape, among a given family, for self-assembling 
the diamond lattice, (ii) Our criterion of 7 /(a*) = 0 is both 
parameter-free and independent of system kinetics, which are 
highly dependent on simulation methods. Nevertheless, we 
find rough agreement between the thermodynamic computa¬ 
tion of the alchemical potential and a measurement of the lower 
critical packing fraction r] c reported in Ref. lt33l . (iii) The fact 
that the optimal particle shape ( a * « 0.37) for diamond as¬ 
sembly at 77 = 0.6 is more tetrahedral than the optimal shape 
for diamond packing (a a = 1/3), but not perfectly tetrahe¬ 
dral {a = 1 ), arises from a competition between two effects. 
Particles must have tetrahedral valence to form the diamond 
lattice, but in the diamond lattice, particles are arranged in an 
alternating motif (Fig. |TJc,d). Shape entropy considerations 
[16] suggest that as the system density is lowered, particles 
must have larger entropic patches fl4l to maintain their emer¬ 
gent valence, as shown in Fig. [2] However, as illustrated in 
Fig.QJl, if the particles are too tetrahedral, then the alternating 
diamond motif leads to overlapping next-to-nearest neighbor 
particles, as shown in Fig. [2] Hence, the optimal truncation 




































6 


of a tetrahedron to self-assemble diamond is more tetrahedral 
than packing would dictate to preserve valence, but not too 
tetrahedral to prevent particles from having alternating valence, 
(iv) We computed the alchemical modulus M a for truncated 
tetrahedra at 77 = 0.6 and a = a*. In future work it would be 
interesting to determine how this modulus varies across system 
density in this system and differs between systems/structures, 
or relates to effects of polydispersity, and how it behaves at 
phase boundaries, (v) The entropic assembly of anisotropic 
hard shapes is driven by emergent valence! 14, 16], manifest¬ 
ing in directional entropic forces (33l . A defining feature of 
emergent behaviors is that their origin is difficult to trace to 
microscopic attributes of the system constituents. Q2| Here, we 
explicitly demonstrate the general principle that it is possible to 
optimize building block attributes, by which we systematically 
control emergent valence, in order to optimally assemble a 
target structure. Moreover, our results suggest a general design 
rule for entropic valence: that as system density decreases, 
entropic patch size fl4l must increase to optimally assemble a 
dense packing phase. This design rule is supported by another 
recent result [43 ] where it was found that for several families 
of dimpled particles, the peak in packing density occurs at an 
entropic patch size that is below the critical size for the onset 
of entropic assembly at low density. This is particularly strong 
evidence for the design rule proposed here because the optimal 
patch size cannot be smaller than the patch size at onset, (vi) In 
practice, the synthesis of anisotropic colloidal particles is often 
driven by a growth process that yields particles in a family 
of shapes. Here we have shown, in an example family, how 
to optimally choose when to terminate that growth process to 
obtain particles for assembling a specific target structure. 


B. Rhombic Dodecahedron 

To (i) understand how to contrast the relative importance of 
different shape modifications of a given shape, and (ii) deter¬ 
mine how this relative importance depends on system density, 
we studied a two-parameter family of truncations of rhombic 
dodecahedra that leave them invariant under the spheric trian¬ 
gle group A 4 ,2,3.[56] The A 4 j 2,3 invariant family of shapes 
is constructed with three families of planes that make up the 
faces of a cube, a rhombic dodecahedron, and an octahedron, 
all oriented to preserve the necessary point group symmetry. 
The rhombic dodecahedron has two different types of vertices: 
four-fold vertices where four planes come together, and three¬ 
fold vertices where three planes come together. Moving the 
planes that make up the faces of the cube towards the origin 
truncates the four-fold vertices, and moving the planes that 
make up the faces of the octahedron truncates the three-fold 
vertices. We performed simulations that examine the effects 
of each type of truncation on a perfect rhombic dodecahedron. 
We parametrize the vertex truncations so that when 014 = 0 
(four-fold vertex truncation) and <^3 = 0 (three-fold vertex 
truncation) the particle is a perfect rhombic dodecahedron. 
Maximal truncation 014 — 1 and <^3 = 0 yields a perfect cube, 
and <X 4 = 0 and <^3 = 1 yields a perfect octahedron. 

We determined how systems of perfect rhombic dodecahedra 


(04 = c ^3 = 0 ) respond to infinitesimal changes in a 3 and 0 L 4 . 
We computed the alchemical potentials /i 4 conjugate to a 4 
(four-fold vertex truncations) and 11 3 conjugate to a 3 (three¬ 
fold vertex truncations) for systems of N = 256 rhombic 
dodecahedra at a series of packing densities 77 between 0.525 
and 0.75 in increments of 0.025 at 014 — <^3 = 0. As shown 
in Fig. [3] we find negative alchemical potentials for both 3- 
fold and 4-fold vertex truncations ( 7 / 3 , /i 4 < 0) at all densities 
studied 0.525 <77 <0.75, implying that both types of vertex 
truncation reduce the free energy of the system. Moreover, we 
find that truncation of the four-fold vertices results in a greater 
reduction in free energy than the three-fold truncation. 

Our computation of the constitutive relations ( 77 ) for rhom¬ 
bic dodecahedra explicitly demonstrates how our methods can 
determine the relative importance of various shape features. 
Determining the most important shape features to control is cru¬ 
cial for anisotropic particle synthesis techniques, and here we 
have demonstrated a general method for solving this problem. 
In addition to providing this general proof-of-principle, our 
results have several specific implications, (i) At all densities 
studied, we observed /i 4 < /i 3 < 0 indicating that both types 
of vertex truncation improve the self-assembly of rhombic do¬ 
decahedra into a face-centered cubic (fee) lattice. Because 
vertex truncation at fixed volume means the particles become 
slightly more spherical, our result suggests that the structure 
is further stabilized by particles exchanging some vibrational 
degrees of freedom for rotational ones. Moreover, (ii) be¬ 
cause /i 4 < /i 3 it suggests that the four-fold vertex truncations 
are more important in restricting the rotational motion than 
the three-fold vertices. There are 8 three-fold vertices and 6 
four-fold vertices in a rhombic dodecahedron, but the centroid- 
to-vertex distance for a four-fold vertex is 4/3 the distance 
for a three-fold vertex. We might suspect that if a vertex type 
sticks out further from the shape, or is greater in number, it will 
provide a greater steric constraint on the microstates available 
to the system. Our result that /i 4 < 773 suggests that for the 
rhombic dodecahedron in an fee lattice, the vertex distance 
is more important than the number of vertices. It would be 
interesting to investigate whether this design rule holds for 
other shapes, or is specific to rhombic dodecahedra. (iii) Be¬ 
cause the slopes of both 77 ( 77 ) curves are positive for 77 > 0 . 6 , 
it suggests that particles give up rotational entropy faster than 
translational entropy as the system density increases. We note 
that the distinction between four-fold and three-fold vertices 
becomes smaller at larger packing fractions, which suggests, 
surprisingly, that as the particles increasingly lose rotational 
entropy the distinction between how they lose it becomes less 
important. It would be interesting to see if this result holds 
more generally in other systems. 


C. Oscillating Pair Potential 

To demonstrate that our alchemy approach is not limited to 
particle shapes, we studied spherical nanoparticles (or point 
particles) interacting isotropically using a truncated, interme¬ 
diate range oscillating pair potential studied in (42), which is 
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Figure 4. Panel a: Alchemical potential {—ft) for the wavenumber k, 
and phase 0 parameters systems of N = 4096 particles interacting 
via a three-well oscillating pair potential over a range of parameters 
that self-assemble an icosahedral quasicrystal. EH The pair potential 
self-assembles icosahedral quasicrystals of three different densities. 
Here, we examine the region of parameter space that self-assembles 
the intermediate density quasicrystal (dashed lines indicate the phase 
boundaries we observed for the self-assembly of the intermediate 
density quasicrystal). Surprisingly, we find that, over the range of 
parameters we studied, in order to thermodynamically improve the as¬ 
sembly of the intermediate density quasicrystal, we are driven toward 
the region of parameter space that is dominated by the self-assembly 
of the high-density quasicrystalline phase. This suggests that the 
optimal choice of parameters to stabilize the intermediate density qua¬ 
sicrystal is buried in a region that will spontaneously self-assemble 
the high density phase instead, and suggests that the intermediate 
density phase will be difficult to stabilize in practice. The insets show 
the bond order diagram and diffraction pattern from a simulation 
snapshot of a 4096 particle system at k = 8, 0 = 0.53. Panel b: 
In the same system we computed the average potential energy per 
particle for different values of k as a function of 0. We see a clear 
decrease in (U) with increasing k and 0. This finding suggests that 
the decrease in free energy with increasing k and 0 shown in panel a 
can be attributed to enthalpic contributions from lower ground state 
energies. 


inspired by Friedel oscillations. It can be written in the form 

U(r) = ^15 + ~3 cos (^( r _ 1-25) — 0) • (11) 

This potential has been recently shown to self-assemble an 
icosahedral quasicrystal for k^T = 0.25 for 0.78 < k < 0.82 
and 0.52 < 0 < 0.55 (42]. The potential is of particular 
interest due to the possibility of realizing it in systems of 
nanoparticles or colloids decorated with appropriate ligands. 
For these computations, we work in units with e = 1. We 
performed simulations of N = 4096 particles using HOOMD- 
Blue.[57|| For full simulation details, see appendix. 

We computed the alchemical potentials conjugate to k 
(wavenumber) and /i^ conjugate to 0 (phase shift) for systems 
of N = 4096 particles interacting via the oscillating pair po¬ 
tential in Eq. fTT) . We studied the pair potential in the range 
of parameter space that was shown previously (42) to self- 
assemble an intermediate density icosahedral quasicrystal. In 
this phase, we find that within the entire parameter range over 
which we were able to reliably nucleate the intermediate den¬ 
sity quasicrystal, both /!& and /i^ are negative. We show this 
explicitly in Fig. [4^ where we form and fi^ into the vector 
ft. We plot — ft, which shows the direction that decreases the 
free energy at a given point in parameter space. 

This result alone does not indicate whether this curious 
behavior is enthalpic or entropic in origin. To understand the 
origin of this decrease in free energy for increasing both k 
and 0, we computed the average potential energy at each state 
point, which is plotted in Fig. B>- We see that at a given k, 
increasing 0 decreases the system’s potential energy, and that 
the potential energy is lower at a given 0 with increasing k, 
which is consistent with the alchemical potential results shown 
in Fig. [4^. This suggests that the effect we observe in Fig. [4^ 
is enthalpic in origin. 

Surprisingly, our result that /i& and /i^ are everywhere nega¬ 
tive suggests that there is not a choice of parameters for which 
ft = 0 (i.e. a local free energy minimum) in the parameter 
regime where the intermediate density quasicrystal is the ther¬ 
modynamically preferred phase. (For an example of a simpler 
case where there is a local free energy minimum in a system 
with isotropic interactions, see appendix.) Rather it suggests 
that, at least for systems of N = 4096 particles, the optimal 
parameter choice for self-assembling the intermediate density 
quasicrystal lies along the boundary separating the assembly 
of the intermediate density quasicrystal and the high density 
quasicrystal, which is the thermodynamically preferred phase 
at higher values of k and 0.(42) 

A general take-away message of Ref. (42 ) is that controlling 
assembly in one-component systems via isotropic interaction 
potentials involves two things. It involves controlling the rela¬ 
tive distances of potential energy minima, which determines 
preferred relative distances between particles. Note that precise 
control over this procedure not straightforward, even at T = 0. 
See appendix for an explicit demonstration in a toy model sys¬ 
tem. However, it also involves controlling the relative depth of 
the minima, which determines the number of particles that sit 
at the preferred relative distances determined by the minima 
locations. Here, we are able to directly compute the effects 














8 



Atomic Matter 

Colloidal Matter 

Anisotropy 

Dimensions 

Proton Number 

Many 

Anisotropy 
Dimension Types 

Discrete 

Discrete, 

Continuous 

Vfilpnpp 

Quantum Mechanics, 


V dlClltC 

Constraints 

Group Theory, 
Fermi Statistics 

Steric 

Number of 
Stable Elements 

-80 

Infinite 


Figure 5. Contrast between constraints on engineering materials with 
atomic elements and colloidal “elements”. Colloidal elements have 
valence that can vary continuously in many different ways. The 
malleability of colloid valence means that constructing a “periodic 
table” for them is inherently difficult. However, we can exploit the 
malleability of colloids to directly probe how particle attributes affect 
structure. 


of changes in potential control parameters on the system free 
energy and we find that they can be detected. In appendix, we 
consider the pattern registration as measured by comparing the 
locations of the potential minima with the radial distribution 
function of the particles, and we find no discernible difference 
across the range of parameters we considered. This suggests 
that our alchemical potential methods are sensitive to system 
behavior that is not easily discernible via conventional analysis. 
We believe this might be of particular value in systems such as 
the oscillating pair potential system where there is a very rich 
bulk phase structure that depends sensitively on the choice of 
potential parameters controlling particle valence (42). 


III. CONCLUDING REMARKS 

We chose families of model systems to demonstrate the 
power of our methods because of their structural complexity 
(the icosahedral quasicrystal), or conceptual complexity (the 
emergent behavior of hard shapes); however our methods can 
be generalized straightforwardly to systems of particles with 
other interactions or shapes, as well as systems with enthalpic 
patches Of) or multiple particle species. Furthermore, though 
our focus was on understanding macroscopic colloidal behav¬ 
ior within a given region of phase space, our methods can 
be applied to the crystallization of other types of matter, e.g. 
polymers, and the study of phase boundaries. One example 
where both are relevant is in the investigation of the polymor¬ 
phism (58] or supramolecular isomerism in crystals of small 
molecules, which is relevant for pharmaceutical applications 

(59). 

Here we focused on solving the problem of determining 
optimal building block attributes for target structures among 
a range of building blocks, from which we were able to ex¬ 
tract design rules for emergent behavior. As a result, most of 
our calculations were of the constitutive relation /i(7V, V, T, a). 
However, for truncated tetrahedra we also considered (Fig.|jJ 
and appendix Movie) the constitutive relation a(N, V, T, fi) 


for particles fixed to sit on a diamond lattice using a simple 
extension of Eq. Q (see appendix for details). All of the 
foregoing discussion concerning interpretation of alchemical 
potentials, including the relation to building block optimality, 
continues to hold, where any quantities computed in extended 
ensembles are, by design, conditional on the externally im¬ 
posed criteria. Using extended ensembles, it is straightforward 
(see appendix for details) to use our techniques for the dis¬ 
covery of building blocks for bulk materials given a suitable 
choice of external design criteria. We leave a full numerical 
investigation of this class of problems to future work. 

Our method for determining optimal building blocks to self 
assemble target structures was based on the desire to make 
quantitative connections between building block attributes and 
bulk behavior. To make our proof-of-principle demonstra¬ 
tion explicit, we ensured that the local minima we identified 
were bona fide global minima by computing exhaustively over 
relevant building block attributes. Rather than compute ex¬ 
haustively as we have here, future investigations should reduce 
computational effort by employing global optimization tech¬ 
niques. Indeed, work aimed at optimizing building blocks 
for bulk attributes has employed genetic or evolutionary algo¬ 
rithms 16QH64) . or gradient descent (65) . Those approaches are 
complementary to the optimization part of the present work in 
three ways: (i) Our approach provides a systematic, rigorous, 
first-principles method for constructing probability distribu¬ 
tions needed to apply the gradient descent method proposed in 
Ref. (65) . (ii) Genetic and evolutionary algorithms are pow¬ 
erful techniques that use external fitness criteria to perform 
non-local optimization. Our approach supplements these non¬ 
local approaches by providing direct, precise measurement 
of the physical response of a system to a local change in the 
attributes of building blocks, (iii) The ability to probe local 
changes in building block attributes is also important because, 
in addition to optimizing attributes, we would like to be able to 
derive generalizable design rules that extend beyond specific 
systems of interest. Here we showed an example of how to 
accomplish this using digital alchemy by showing that dense 
packing arguments for anisotropic shapes can be extended to 
lower density by increasing the size of entropic patches. We 
believe that a combination of the methods we present here with 
existing techniques (6014631 [65] will provide a powerful tool 
set for materials design. 

Finally, our digital alchemy method shows how to phrase a 
generic class of relationships between building block attributes 
and bulk behavior for colloidal materials. For colloids, the fact 
that valence can vary in many ways (sometimes continuously) 
along several different anisotropy dimensions iE) means 
that it is not possible, even in principle, for a periodic table of 
colloidal elements to be as succinct as the atomic periodic table 
(see Fig. [5]). However, like the atomic periodic table relates 
atomic valence to bulk behavior, we have shown that it is possi¬ 
ble to relate colloid valence to bulk behavior. Indeed, because 
colloid valence is so malleable, we have shown that building 
block property-bulk behavior relationships for colloids can be 
quantitative in a way that is not possible for atoms. In effect, 
whereas quantum mechanics dictates that the atomic periodic 
table is complete and succinct, but heuristic, the outlines for a 
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periodic table of colloidal elements suggested by this work are 
that it is complex and many-dimensional, but also quantitative, 
and richly predictive. 
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Appendix A: Supplementary Theory and Methods 


only. Defining 

e -pF s J[ dq ][ d Q] e -0u {ai} ( q ,Q) } (A3) 

we have, up to irrelevant constants, the thermodynamic poten¬ 
tial 

N 

PF = log det (/{„,}) + /3F . (A4) 

We compute alchemical potentials by differentiating this ex¬ 
pression with respect to the a^. Even if the particle mass and 
volume are fixed, the first term depends on the particle shape, 
and so we need to compute the moment of inertia tensor of our 
particles. 

This term does not contribute to the computation of isotropic 
(spherical particles), however it is important for polyhedral par¬ 
ticles. Our computation used four steps, (i) We compute the 
moments of inertia by identifying all of the faces of the polyhe¬ 
dron. (ii) We do a fan decomposition of the faces into triangles, 
(iii) We use the point at the origin with the triangulation of 
each of the faces to decompose the polyhedron into a set of 
tetrahedra. (iv) We use standard formulae to compute the in¬ 
ertia tensor of the tetrahedron. ll66l We checked that our code 
was correct by using it to compute moments of inertia for 
several known shapes. As an additional check, all shapes we 
considered are invariant under triangle group symmetries. Via 
Schur’s lemma (see, e.g. l67l ). their moments of inertia tensors 
must be proportional to the identity matrix, and we checked 
that our code gave results consistent with this up to machine 
precision. 


1. Moments of Inertia 


2. Numerical Evaluation of Alchemical Potential 


Our calculation of the alchemical potential that determines 
the constitutive relation Pi({oa}) depends on the moment of 
inertia tensor of the system. To see this, we write the partition 
function for the ensemble with fixed as 

Z({di}) oc J [dp\[dL][dq][dQ}e~ 0H , (Al) 

where q are particle positions, Q are particle orientations, p 
are conjugate momenta, L are angular momenta, we have 
suppressed particle indices, and for simplicity we are working 
in an ensemble with fixed volume and number of particles. 
The following discussion is straightforward to extend to other 
ensembles. 

Starting with Eq. ( |A1| ), we integrate over the momenta and 
angular momenta, which gives (in three spatial dimensions) 

Z({di}) oc/? _3iV TO 3JV/2 (det(/ {a . } )) iV/2 

r , s (A2) 

/ [dq][dQ]e- pu ^ q ^ . 

For simplicity, we will concern ourselves with changes in 
alchemical parameters that leave the particle mass and volume 
invariant so that we are isolating the effects of changes in shape 


To evaluate the contribution of the configuration integral, 
F, to the alchemical potential we use a variant of the Bennett 
acceptance ratio method. |48) 

We evaluate the expression 


1 dF 
“ Ndal 


(A5) 


using finite differences. For brevity, we will give formulae for 
a single a\ the extension to multiple a is straightforward. For 
some finite h, we estimate 


fa ~ iT2^ F ^ a+vh) > < A6) 

V 


where, and v are appropriate constants for some finite dif¬ 
ferencing scheme 1688 . To suppress numerical errors we used 
a symmetric four-point scheme for calculations involving trun¬ 
cated tetrahedra, the oscillating pair potential, and the 2D 
Lennard-Jones-Gauss system (see below); for rhombic dodeca- 
hedra we used a one-sided four point scheme. 

For each state point a , we independently equilibrated several 
copies of the target crystal lattice at each nearby ol v = a + vh. 
For the N = 216 (h = 4 x 10 -4 ) systems of truncated 
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tetrahedra, the N = 4096 oscillating pair potential systems 
(h k = 3 x 10“ 5 , h ^ = 10“ 5 or 2 x 1(T 5 ), and the N = 1024 
(h = 10 -5 ) systems of Lennard-Jones-Gauss particles we 
self-assembled crystals from the fluid; for the N = 1000 
(h = 10 -4 ) systems of truncated tetrahedra and the N = 256 
systems of truncated rhombic dodecahedra (h = 2 x 10 -3 ) 
we constructed the crystal directly. Without loss of general¬ 
ity, we labelled a particular v as v§. For the hard truncated 
tetrahedra and rhombic dodecahedron systems, for all v ^ z/ 0 
we repeatedly sampled states from the equilibrium distribution 
using standard Metropolis Monte Carlo techniques (see, e.g. 
(54) k and computed the probability, according to the Metropo¬ 
lis criterion (69), of accepting a trial Monte Carlo move of 
the state from the ensemble with v to the ensemble with z/q- 
Similarly, we repeatedly sampled states from the equilibrium 
distribution of the system at z^o, and computed the probability, 
according to the Metropolis criterion, of making a trial Monte 
Carlo move of the state from ensemble vo to each of ensembles 
at the other v. For the oscillating pair potential and Lennard- 
Jones-Gauss systems, we repeatedly sampled configurations 
from equilibrium NVT molecular dynamics trajectories cou¬ 
pled to a Langevin thermostat in HOOMD-Blue. (57) At each 
sample point we computed the potential energy in the system 
at a given v as well as what the potential energy would be if 
the particles maintained all of their positions but interacted 
with a potential z/. By recording these potential differences 
we constructed the probability, according to the Metropolis 
criterion, of making a Monte Carlo move from one value of a 
to another. 

To show how this computes the alchemical potential, we note 
that combining Eq. ( |A5| ) with Eq. (A6| ) and exponentiating we 
have 

e -hN/3fi ^ e ~ a„(3F(a+vh) _ TT e —+vPF(a+vh) . (Nl) 


Decomposing the thermodynamic potential into the kinetic and 
configuration components using Eq. (A4| ) we have 


e ~hN (3 [a 


e 71°S det I (a+vh)—j3F (a+vh)] 
v 


(A8) 


We then use detailed balance to write the configurational part 
of the free energy F at a + vh in terms of the value at a + v 0 h 
to get 

e ~hNf3n ^ logdet I (a+vh)—(3 F(a+v 0 h)] 

v 

/ p(a + i/oh\a + vK) 

\p(a + vh\oi + z/q h) 

where p are the relevant transition probabilities we compute 
with the Bennett acceptance ratio method. We note that, since 
we are evaluating a first derivative of F, Yh v lv = 0 so that 

-hNPn ~ TT„7^ flog det I (a+vh) ( + 1++ + <+) 

\P+ + vh\a + u 0 h) J 

(A10) 




Figure 6. The functional form of the Lennard-Jones-Gauss potential 
(Eq. (|CTJ) for = 3/4, cr G /a L j = a/2/10, and r 0 /cr L j in the 

range for a two dimensional system to form a square lattice. 


Finally, taking the logarithm of both sizes we get that 

~r F, lo g det/ (« + vh ) 

\ ( , ,, , h \ (All) 

1 ^ P\ a + z^o^|c^ + vh) 

Nh “ p(a + vh\oi + z / 0 h) 

Note that we did not attempt to determine general criteria for 
choosing optimal finite differencing schemes, i.e. the value of 
h or the order of the method. Two considerations that arise are 
that if h is too large, the finite differencing error is large, and a 
higher order method must be used. However, if h is too small, 
this affects the autocorrelation time of the transition probabil¬ 
ities which, necessitates more extensive sampling. Figuring 
out how to do this optimally for this type of computation is an 
interesting problem for future work. 


3. Error Estimation 


Numerical evaluation of the expression Eq. ( |A11| ) involves 
both systematic and statistical error. Statistical error comes 
from the estimation of the transition probabilities, and is given 
by 


(5(/3/i)stat — 


1 \ /8 p{oL -f- V{)h\ cm, H - vh) 

Nh V p(a + voh\a + vh) 


Sp(a + vh\a + vq h) \ 
p(a + vh\a + vq h) ) 


(A 12) 


Systematic error can arise from the calculation of the mo¬ 
ment of inertia tensor. In practice, the moment of inertia tensor 
can be computed to machine precision, so this contribution 
is negligible. Further systematic error can arise if systems 
are equilibrated at slightly different densities. To see how 
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Figure 8. We examine pattern registration of self-assembled square 
lattices in 2D as shown through the radial distribution function (g(r), 
darker curve) with the Lennard-Jones-Gauss pair potential (lighter 
curve, minima indicated with darker marks) that induces it to assem¬ 
ble, for various values of the parameter a = ro /ctl j (see Eq. in 
SI text). In these plots ro is less than the critical value of a* where 
the alchemical potential fi vanishes. 


Figure 7. Constitutive relation fi(a) (■, top plot), where a = to/ctl j 
describes the potential well position, for a two-dimensional Lennard- 
Jones-Gauss system of 1024 particles at k^T — O.I^lj (top). In 
this range of a the system forms a square lattice. The alchemical 
potential vanishes at a* ~ 1.38, indicating that that is the optimal 
well-location for self-assembling a square lattice. This accords with a 
calculation showing that the potential energy is minimized at this well 
position both k^T — 0.1 (♦, middle plot) and at k^T = 0 (solid line, 
middle plot), but is lower than the naive geometric estimate of that 
comes from pattern registration of the potential minima to the ratio of 
nearest-neighbor to next to nearest-neighbor distances in the square 
lattice (>/2, intercept in lower plot). 

this arises, consider the differential of the free energy for our 
systems, which has the form 

dF = iiNda - PdV . (A13) 

where N is the number of particles, and for simplicity we 
are considering a single alchemical parameter a and working 
in ensembles with fixed volume. To compute the alchemical 


potential fi we need to differentiate the free energy at fixed vol¬ 
ume. However, there can be a variation in the system volume, 
which in our systems is reflected in a change in the packing 
fraction rj, defined by 



T] 


where £ 3 is the volume of a particle. Differentiating gives 

N£ 3 

dV = - ^dg . (A15) 

7 / 

So the differential for the free energy is given by 

NPi 3 

dF = fiNda H- A — dr] . (A16) 

T] Z 

That means that to get a measurement of the alchemical poten¬ 
tial, we need to have 

P£ 3 

fida >> — dr] (A17) 

r] z 
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Figure 9. We examine pattern registration of self-assembled square 
lattices in 2D as shown through the radial distribution function (g(r), 
darker curve) with the Lennard-Jones-Gauss pair potential (lighter 
curve, minima indicated with darker marks) that induces it to assem¬ 
ble, for various values of the parameter a = ro / cr l j (see Eq. <CTJ in 
SI text). In these plots ro is greater than the critical value of a* where 
the alchemical potential fi vanishes. 


systems. We tested controlling for this error in for four-fold 
vertex truncations at a density of 77 = 0.625 where we observed 
since we observed the maximal sensitivity to packing fraction 
for that type of truncation at that density. We controlled for 
the error by independently equilibrating systems at packing 
fractions r\i distributed near the desired packing fraction rj. We 
obtained a large number of independent samples of the tran¬ 
sition probabilities for the various 77 *. Given these transition 
probabilities, we used regression (weighted by the errors in 
each of the transition probabilities at rji) to estimate the value 
of p at 77 . We used the error from the regression estimate as the 
statistical error in /3/x in Eq. ( |A12| ). To well-within the error 
bounds we found no difference between the results obtained 
with densities randomly spread about the desired value, and 
those set to the desired value to machine precision. 


4. Numerical Limits on Determining Optimal Building Blocks 


We note that in order to determine the roots of the consti¬ 
tutive relation fi{a) that determine optimal particle configura¬ 
tions to arbitrary accuracy runs up against the limit Eq. ( |A17| ). 
In practice we note that reasonable computations on modern 
hardware allow us to bound optimal shapes, for example, to 
differences in morphology that are imperceptible. 


to safely control the error arising from changes in packing 
fraction. We did this in two ways. For the systems of N = 
1000 truncated tetrahedra and N = 256 rhombic dodecahedra 
we built perfect crystals at the desired packing fraction and 
then thermally equilibrated them. E.g. for the systems of N = 
256 rhombic dodecahedra, this fixes the variation in packing 
fraction to < 10 -6 , and so the right hand side of Eq. ( |A17| ) is 
~ 10 -5 . On the left hand side da = 2 x 10 -3 , and /i ~ 10 _1 , 
which indicates that the spread in packing fraction contributes 
a systematic error on the order of 10 %, which does not affect 
any of our conclusions. Systems of N = 1000 tetrahedra are 
similar, but in that case we also find agreement between our 
alchemical potential computation and the computation with 
fluctuating shape (c.f Fig. 2d,e, main text). For the systems 
of TV = 216 tetrahedra, we self-assembled crystals to within 
1 % of the desired packing fraction and then controlled for the 
errors in packing fraction statistically, and again the results of 
that computation agree with the computation with fluctuating 
shape. 

For further confirmation, we note that our computation in¬ 
volves making trial moves between different types of particles 
in the same system. That means that when we compute the 
probability of going from a a + 8 uh, the volume of the 
system is invariant. And, similarly, when we compute the prob¬ 
ability of going from a+5ish a , the volume is also invariant. 
However, there is a small difference in density between the two 


5. Reparametrization and the Alchemical Constitutive 
Relation 

Our calculations of optimal shapes and interactions for 
target structures were carried out by choosing a particular 
parametrization of the particle shape or interaction potential. 
This choice was not unique. For simplicity, consider a sin¬ 
gle parameter family of shapes or interaction potentials, and 
reparametrize them by a = a(a'). Under this reparametriza¬ 
tion, the constitutive relation becomes 


y!{ot) = 7 ^(a) • (A18) 

From this form we note that if the reparametrization is mono¬ 
tonic da'/da > 0 , then sgn (//(a/)) = sgn (/i(a)), which 
means that reparametrization will not change the direction 
of the change a that corresponds to a decrease in free en¬ 
ergy. Moreover, if da'/da % 0, any roots a* of the original 
constitutive relation /i{a) = 0 will coincide with roots of 
the reparametrized constitutive relation ti'(a') according to 
a'* = a'(a*), which means that particle shapes, e.g., deter¬ 
mined to be optimal from our alchemical potential calculations 
are optimal regardless of the way in which particle shape is 
parametrized. 
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Appendix B: Field Directed Alchemy: Design 
1. Extended Ensembles 

In the main text, we are concerned with optimizing among 
known building blocks for a target structure, so we began with 
Eq. (2) (main text) to derive and interpret unbiased statistical 
ensembles. In order to design building blocks for a structure 
for which we do not have a set of a priori candidates, it is 
necessary to modify Eq. (2) (main text) so that the statistical 
ensembles are biased to form a structure with a desired property. 
To make this explicit, we suppose there is some quantity A that 
when evaluated on the design structure takes the value (A). 
The ensemble for this system can be found by maximizing the 
entropy 


up to irrelevant overall multiplicative constants. Detailed bal¬ 
ance requires that for a Markov chain Monte Carlo integration 
to converge to Eq. ( |B4| ) the ratio of the probability of making a 
move from a shape aq to a shape 0^2 IIi ^2 to the probability 
of the reverse move n 2 _^i is equal to the ratio of probabilities 
of being in those states 


n 2 ^i _ 7Ti 

III —>-2 7 T 2 


Using Eq. ( |B4| ) we have 


n 2 ^i _ (det(/ ai )) JV/2 

“ (det (I aa )) N/2 e -/3(^2-MJV«2—AA) ' 


(B5) 


(B 6 ) 


s = - \*<J M 7 ^) - /? (ncrH - (E) 

<J 

-^2^iN(Tr a ai - {on )) 
i 

- A(tt < t A - (A))) , 

(Bl) 


with respect to 7 r a , where A is a Lagrange multiplier. This 
yields the extended partition function 






a 


(B2) 


We performed simulations of truncated tetrahedra with fluctu¬ 
ating shape at zero alchemical potential fi = 0 , in an external 
field A that forces the particles to sit in an Einstein crystal 
with spring constant k (measured in units of k^T/£ 2 ). Al¬ 
chemical moves were performed at fixed particle position and 
orientation, which allows us to simplify Eq. ( |B 6 | ) to 


n 2 ^i _ (det {Ion)) ^ -P(u ai -u a2 ) 
ni ^2 " (det (I a2 )) N/2 


(B7) 


To satisfy detailed balance, we take the Metropolis l69l crite¬ 
rion as 


Note that A, which is a Lagrange multiplier, determines the 
strength of the coupling to the external field. We note that if 
A is positive (negative), the system is driven toward particle 
configurations that favor increasing (decreasing) A, which 
allows one to design both toward or away from a structural 
characteristic encoded in A. Moreover, one could certainly use 
multiple structural characteristics with the aim of arriving at 
building blocks suitable for reconfigurable structures, or the 
suppression of some particular polymorph. 

2. Fluctuating Shape: Detailed Balance and Interpretation 

Here we show how to satisfy detailed balance in systems 
with fluctuating shape aimed at directly solving the problem 
of particle design. We begin with the generalized partition 
function from Eq. ( |B2| ) 

Z = jda[dp}[dL][dq\[dQ}e-^ H -^ Na - x ^ , (B3) 

where we have now explicitly included the measures for the 
integration over the particle momenta p , angular momenta L , 
positions q, and orientations Q , and for notational simplicity 
we use a single anisotropy parameter (the generalized form is 
straightforward). Taking the Hamiltonian from Eq. (1) (main 
text), we perform the quadratic integrals over p and L to get 

Zoc Jda[dq][dQ] (det(/„)) JV/2 e _/3 ^ c/ “ _ ' lJVa_AA ^ , (B4) 


IW = min (l, (B8) 

( (det(I a2 )) N/2 J 

In Fig. le (main text) we report a(p = 0) for truncated 
tetrahedra at a packing density of 77 = 0.6 using an externally 
imposed field that puts the particles in an Einstein diamond 
crystal with spring constant k. We show that a(p = 0) in¬ 
creases with k (i.e. particles become more tetrahedral for large 
k). This result has two implications, (i) It shows that we can 
optimize particle shape not only for the diamond structure, but 
that we can optimize particle shape for a diamond structure 
at a fixed density with a stiffness that is determined by the 
stiffness of the Einstein crystal we impose externally. This sug¬ 
gests, more generally, that digital alchemy can optimize both 
structures and properties of structures, (ii) It shows that one 
effect of making particles more tetrahedral is that they form 
a diamond lattice that is more stiff at fixed density. At large 
spring constants we observed a « 0.5, however, we failed 
to observe the spontaneous assembly of diamond lattices in 
our simulations at truncations that were this small. Moreover, 
our alchemical potential calculations show that as the vertex 
truncation of tetrahedra decreases past the optimal value, par¬ 
ticles lose entropy in the diamond lattice. These two results 
together show that though it is possible to tune the stiffness 
of the crystal, which might be surprising because entropy is 
the only governing property in these systems, the range over 
which additional properties can be tuned while still ensuring 
spontaneous self-assembly is limited by kinetic factors. 
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Figure 10. For the oscillating pair potential at fixed k (Panel a: 
8.0, b: 7.8, c: 8.2) we show the pattern registration as measured by 
coincidence of peaks in the radial distribution function for a snapshot 
of the icosahedral quasicrystal (g(r), darker curve), with the pair 
interaction potential (U (r, k, 0), lighter curve). We see no discernible 
difference in pattern registration over this range of parameters. 

Appendix C: Toy Model: Lennard-Jones-Gauss Potentials in 2D 

As a consistency check, and as a non-trivial check on our 
analysis routines, we apply them to a system for which we 
can determine optimal microscopic parameters for the a given 
macroscopic state through direct calculations. 

We study the effects of a one-parameter deformation of the 
relative position of potential minima in 2D Lennar d-Jones - 
Gauss systems [24) via molecular dynamics simulations with 
HOOMD-Blue 11571 . The potential is given by 

Vug = £u ((^)“ - 2 (^)“) - 

(Cl) 



Figure 11. The variation in the oscillating pair interaction potential 
with <j) at fixed k — 8.0 is small, yet leads to substantial gradients 
in the free energy (Fig. 4a, main text), and detectable differences in 
average potential energy per particle (Fig. 4b, main text), that are not 
discernible from the pattern registration depicted in Fig.fTOfr. 

The phase diagram of this system has been previously deter¬ 
mined in [24 , 25]. We use potential parameters as in prior 
published work l25l : £q = |^lj and cfq = lj, and we 
work at T = £lj/10. Finally, we define a = ro/oTj, and 
take it to be in the appropriate range to self-assemble a square 
lattice. 

We studied the constitutive relation between (i and a = 
^o/°lj for the Lennard-Jones-Gauss system in 2D. In Fig. [7] 
we plot the alchemical potential for a system of N = 1024 
Lennard-Jones-Gauss particles at T = O.I^lj for 1.28125 < 
a < 1.46875 where we observed the formation of a square 
lattice. To determine the root of the constitutive equation, we 
performed a linear fit 

— oc (a - a*) , (C2) 

£lj 

where we find a* = 1.383 ±0.001. For a square lattice, the 
ratio of the distance between first neighbours and second neigh¬ 
bours is y/2. We numerically computed the locations of the 
minima of the Lennard-Jones-Gauss potential for ro in the 
range where we observed formation of the square lattice, and 
found that ratio of the second minimum r* ; to the first mini¬ 
mum r* 1 ) was y/2 when a ~ 1.4267, which is well above the 
optimal value indicated by the alchemical potential calculation. 
Because the temperature is low, we expect the free energy to be 
dominated by the potential energy, so we computed the average 
potential energy, and found that there was a potential energy 
minimum between 1.359375 < a < 1.390625 in agreement 
with our alchemical potential calculations. 

Our alchemical potential calculation for the oscillating pair 
potential showed that it is possible to detect effects that are not 
readily apparent by examining the pattern registration between 
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g(r) and the pair potential, as shown in Fig. 10 However, a sur¬ 
prising result of the computations for the Lennard-Jones-Gauss 
system in 2D is that even in simple cases where pattern regis¬ 
tration effects are discernible, the optimal pattern registration 
is not what one would anticipate from naive guessing. 

In the Lennard-Jones-Gauss system, computations were per¬ 
formed at relatively low temperatures, T = O.I^lj, and we 
would expect that at sufficiently low temperatures, the free 
energy of the system is dominated by the potential energy. In 
this case, we expect that the optimal a * we’ve determined 
above coincides with the value of a with the lowest ground 
state energy. For the range of a that self-assemble the square 
lattice, we compute the ground state energy as a function of 
both a and the lattice spacing. For each a we found the lattice 
spacing with the lowest energy to get the ground state energy 
of the square lattice as a function of a. From this curve, we 
found the value of a with the minimum ground state energy to 
be a = 1.38342, which accords very well with the a* we com¬ 
puted using the alchemical technique at T = O.I^lj (c.f Fig. 
[7] top panel and middle panel). However, we note that Figs. [5] 
and [9] show the surprising result that neither result corresponds 
with the naive ansatz of a a = 1.4267 which that comes from 
fixing the relative distance between the first and second min¬ 
ima of the Lennard-Jones-Gauss potential /r^ = y/2, 
which we might expect to optimal because it coincides with 
the appropriate distances for the square lattice in 2D (see Fig. 
[7J lower panel). It is possible that this discordance between 
the naive ansatz from pattern registration considerations, and 
the thermodynamically optimal pair potential might differ in 
other systems, which could be an important consideration in 
DNA-mediated nano-particle superlattice assembly. f70U80l 


Appendix D: Oscillating Pair Potential 
1. Simulation Protocol 

To evaluate the alchemical potentials conjugate to k and (j> 
for the oscillating pair potential system J42], we performed MD 
simulations of N = 4069 particles using HOOMD-Blue lf57lL 
using a tabulated potential to directly reproduce the simulation 
technique employed in l42l . The system size of N = 4096 
particles was chosen so that sufficiently large changes in the 
k and <fi parameters of the potential could be made, and still 
replicate the simulation protocol followed in [[42]]. First, at 
each state point (fc, 0) we performed NVT simulations with a 
cooling schedule that was linear in temperature from an initial 
temperature of3eto0.525e over 5 x 10 7 time steps, and then 
further from 0.525e to 0.25e over 5 x 10 7 more time steps, to 
reach a supercooled fluid. From the supercooled fluid snap¬ 
shots, we launched several NVT simulations to nucleate the 
quasicrystal, in each case seeding the random number gener¬ 
ator of the Brownian integrator with a different integer. As 


a consistency check we determined that on the boundaries of 
the stable range for the intermediate density quasicrystal we 
observed a substantial fraction of events in which we observed 
the nucleation of structures consistent with low density or high 
density quasicrystal where appropriate, which suggests that 
due to Lyapunov instabilities, our procedure leads to uncorre¬ 
lated bulk structures over the whole range. We allowed each 
simulation to run for up to 3 x 10 8 time steps, checking the 
potential energy every 10 6 time steps. Based on empirical crite¬ 
rion of U/Ne < 0.25 we determined that the quasicrystal had 
nucleated or was about to nucleate, and then ran the simulation 
for a further 5 x 10 6 time steps. For each putative nucleated 
quasicrystal, we examined the structure and compared it with 
the structures reported in [42j. Note that, likely due to the 
relatively small number of particles (N = 4096) we found the 
region of self-assembly of the intermediate density quasicrystal 
was smaller that that reported based on simulations of larger 
systems in ll42l . For each intermediate density quasicrystal, we 
equilibrated for 1.3 x 10 8 time steps, and then over a period 
of 1.2 x 10 8 time steps we stored a snapshot of the system 
every 4 x 10 6 time steps. All of the above time scales were 
determined to ensure decorrelation based on autocorrelation 
measurements of the potential energy, and observing diffusion 
of the quasicrystal in the simulation box. From these many 
independent samples of the structure at a given (&,</>), we ran¬ 
domly selected snapshots to re-equilibrate at nearby (fc',0) 
or (fc, (j)') according the finite differencing scheme described 
above, which we did over 10 7 time steps, again chosen to 
ensure statistical independence based on measurement of po¬ 
tential energy correlation. We then repeatedly sampled the 
potential energy every 10 5 time steps over a further 10 7 simu¬ 
lation time steps in order to estimate the probability of making 
ghost Monte-Carlo moves between different values of k and 0 
as described in detail above. At each nearby value of (fc, (j>) we 
obtained several independent estimates of the transition proba¬ 
bilities, and estimated the error from the standard deviation of 
the distribution of the independent estimates. 


2. Pattern Registration in the Oscillating Pair Potential 


We computed the radial distribution function g(r) for equi¬ 
librated snapshots of the oscillating pair potential system at 
k^T = 0.25e and various values of k and and compared it 
with the pair interaction potential U (r) in Fig. 10 to determine 
whether the decrease in free energy we found in the alchemical 
potential calculation (Fig. 4a, main text) could be detected in 


the pattern registration. Fig. [10]shows no clearly discernible 
difference in the pattern registration. Furthermore, in Fig. 03 
we plot potential energy difference as a function of r at fixed 
k = 8.0 between 0 = 0.55 and </> = 0.52, and see that over 
the range of the oscillating pair potential, the differences are 
less than 3% of e. 
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